22.46 Procesamiento Adaptativo de Señales Aleatorias

Trabajo Práctico N° 2

Filtrado lineal óptimo

En este trabajo se modela la respuesta al impulso $h(n)$ de un sistema parlante-habitación-micrófono, asumiendo régimen estacionario. Para lo cual se analizan diversas alternativas de señales de sonido con las cuales excitar al sistema para medir el comportamiento del mismo a la salida. Se emplea esta información en un esquema de identificación de sistemas para modelar el comportamiento de este sistema físico desconocido. Finalmente, se busca aplicar el modelo obtenido para recuperar la señal.

Grupo N° 1

• Davidov, Gonzalo Joaquín
• Farall, Facundo David
• Kammann, Lucas Agustín
• Trozzo, Rafael Nicolás

1. Selección de señales de excitación $x_j(n)$

En esta sección se generan diferentes señales de sonido $x_j(n)$ que serán utilizadas para excitar al sistema parlante-habitación-micrófono, con el objetivo de medir la respuesta $y_j(n)$ de dicho sistema. Para todos los casos, se utilizará una frecuencia de muestreo $f_s = 48kHz$. Las muestras de sonido tendrán un intervalo de 10 segundos.

Por otro lado, se busca verificar visualmente que las señales de entrada y salida muestreadas se encuentren sincronizadas o, al menos, sin infringir la causalidad del evento.

import matplotlib.pyplot as plt

from scipy import signal

from scipy.io.wavfile import write

from scipy.io.wavfile import read

import numpy as np

import IPython

import seaborn as sns

sns.set_style('dark')
sns.set_context('paper')
sns.set_palette('flare')

def plot_input_and_output(x, y, size, title):
    """ Plot the input and the output of the physical system to visually compare them
        @param x Input samples
        @param y Output samples
        @param size Amount of samples to show
        @param title Title of the plot
    """
    # Create the subplots to layout the plots
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(20, 8))
    fig.suptitle(title, fontsize=20)
    
    # Plot the input samples
    ax1.plot(x[:size], color='blue')
    ax1.set_ylabel('$x(n)$', fontsize=18)
    ax1.set_xlabel('$n$', fontsize=18)
    ax1.grid()
    
    # Plot the output samples
    ax2.plot(y[:size], color='red')
    ax2.set_ylabel('$y(n)$', fontsize=18)
    ax2.set_xlabel('$n$', fontsize=18)
    ax2.grid()

# Define the time interval
fs = 48e3
time_interval = 10
t = np.linspace(0, time_interval, int(time_interval * fs))

1.1. Señal de voz

Para modelar el sistema físico desconocido (en este caso la habitación), es necesario que exista correlación entre la señal de entrada y la señal deseada, y sólo se podrá ajustar en aquellas frecuencias presentes simultáneamente en ambas señales. Por ende, en primer lugar es necesario que las señales utilizadas tengan el mayor contenido espectral posible. No se podrá estimar el comportamiento del sistema para aquellas frecuencias en donde no hay contenido espectral en la señal de prueba usada. En conclusión, es necesario una señal con la mayor cantidad de frecuencias posible, idealmente todas (las audibles). De esta forma, se podrá excitar al sistema en todos los modos posibles.

Naturalmente, la respuesta impulsiva se podría medir a partir de excitar al sistema con un impulso o delta de dirac, y justamente presenta un espectro constante o plano, con contenido en todas las frecuencias. Por ende, la señal a escoger podría además tener características impulsivas. Esto es, aparición de tonos repentinos, que permitirían al estimador obtener información del sistema en los modos excitados por tales tonos. Finalmente, otro aspecto de vital importancia para el experimento, es garantizar que los modos del sistema se vean excitados con una cantidad de energía apropiada. Así, la estimación en función de lo percibido en la salida será mejor.

1.1.1. Señal entrante $x_1(n)$

fs_x_voice, x_voice = read('../resources/voice.wav')

IPython.display.Audio(x_voice, rate=fs_x_voice)

1.1.2. Señal saliente $y_1(n)$

fs_y_voice, y_voice = read('../resources/Voice-quiet.wav')

y_voice = y_voice[:len(x_voice)]

IPython.display.Audio(y_voice, rate=fs_y_voice)

1.1.3. Comparación $x_1(n)$ vs $y_1(n)$

plot_input_and_output(x_voice, y_voice[25000:], 10000, 'Señal de voz')

1.2. Señal de música

Para la elección de la señal de música, se optó por utilizar la canción Bohemian Rhapsody de Queen. La razón por la cual se escogió esta canción en el entorno desde el minuto 3:45 al minuto 3:55, es porque en esa región se considera que se cumplen los criterios establecidos anteriormente (en la sección 1.1. Señal de voz).

1.2.1. Señal entrante $x_2(n)$

fs_x_music, x_music = read('../resources/music.wav')

IPython.display.Audio(x_music, rate=fs_x_music)

1.2.2. Señal entrante $y_2(n)$

fs_y_music, y_music = read('../resources/Music-quiet.wav')

y_music = y_music[:len(x_music)]

IPython.display.Audio(y_music, rate=fs_y_music)

1.2.3. Comparación $x_2(n)$ vs $y_2(n)$

plot_input_and_output(x_music, y_music, 10000, 'Señal de música')

1.3. Señal rectangular

1.3.1. Señal entrante $x_3(n)$

f = 100
x_square = signal.square(2 * np.pi * f * t)

write('../resources/square.wav', int(fs), x_square.astype(np.float32))

IPython.display.Audio(x_square, rate=fs)

1.3.2. Señal saliente $y_3(n)$

fs_y_square, y_square = read('../resources/Square-quiet.wav')

y_square = y_square[:len(x_square)]

IPython.display.Audio(y_square, rate=fs_y_square)

1.3.3. Comparación $x_3(n)$ vs $y_3(n)$

plot_input_and_output(x_square, y_square, 10000, 'Señal rectangular')

1.4. Barrido lineal de frecuencia

1.4.1. Señal entrante $x_4(n)$

f0 = 20
f1 = 20e3
x_linear_sweep = np.sin(2 * np.pi * f0 * t + np.pi * ((f1 - f0) / time_interval) * t**2)

write('../resources/x_linear_sweep.wav', int(fs), x_linear_sweep.astype(np.float32))

IPython.display.Audio(x_linear_sweep, rate=fs)

1.4.2. Señal saliente $y_4(n)$

fs_y_linear_sweep, y_linear_sweep = read('../resources/Linear-sweep-quiet.wav')

y_linear_sweep = y_linear_sweep.astype(float)
# Normalization
y_linear_sweep /= y_linear_sweep.max()

IPython.display.Audio(y_linear_sweep, rate=fs_y_linear_sweep)

1.4.3. Comparación $x_4(n)$ vs $y_4(n)$

plot_input_and_output(x_linear_sweep, y_linear_sweep[:], 100000, 'Barrido lineal de frecuencia')

1.5. Barrido exponencial de frecuencia

1.5.1. Señal entrante $x_5(n)$

f0 = 20
f1 = 20e3
k = (f1 / f0)**(1 / time_interval)
x_exponential_sweep = np.sin(2 * np.pi * f0 * (k ** t - 1) / np.log(k))

write('../resources/exponential_sweep.wav', int(fs), x_exponential_sweep.astype(np.float32))

IPython.display.Audio(x_exponential_sweep, rate=fs)

1.5.2. Señal saliente $y_5(n)$

fs_y_exponential_sweep, y_exponential_sweep = read('../resources/Exp-sweep-quiet.wav')

IPython.display.Audio(y_exponential_sweep, rate=fs_y_exponential_sweep)

1.5.3. Comparación $x_5(n)$ vs $y_5(n)$

plot_input_and_output(x_exponential_sweep, y_exponential_sweep[50000:], 100000, 'Barrido exponencial de frecuencia')

1.6. Ruido blanco gaussiano

1.6.1. Señal entrante $x_6(n)$

mu = 0
sigma = 1
x_noise = np.random.normal(loc=mu, scale=sigma, size=len(t))

# write('../resources/white_noise.wav', int(fs), x_noise.astype(np.float32))
fs_x_noise, x_noise = read('../resources/white_noise.wav')

IPython.display.Audio(x_noise, rate=fs)

1.6.2. Señal saliente $y_6(n)$

fs_y_noise, y_noise = read('../resources/Noise-quiet.wav')

y_noise = y_noise[:len(x_noise)]

IPython.display.Audio(y_noise, rate=fs_y_noise)

1.6.3. Comparación $x_6(n)$ vs $y_6(n)$

plot_input_and_output(x_noise, y_noise, 10000, 'Ruido blanco gaussiano')

2. Análisis de las excitaciones $x_i(n)$ y las respuestas $y_j(n)$

En esta sección, realizamos un análisis espectral de las señales de excitación y las respuestas grabadas del sistema. Inicialmente, se buscó graficar los periodogramas de las señales de excitación, para estudiar la conveniencia de cada señal. Se agregaron, luego, los periodogramas de las señales de salida y espectrogramas, para tener un análisis más completo y útil para los ejercicios posteriores.

from scipy import fft

def plot_periodograms(x, y, fs, title):
    """ Plot the power spectrum of the input signal and the desired response signal
        @param x Input signal
        @param y Desired response signal
        @param fs Sampling frequency
        @param title Title of the plot
    """
    # Estimate the power spectrum density
    fxx, Rxx = signal.welch(x, fs, nperseg=4096, scaling='spectrum', window='hanning')
    fyy, Ryy = signal.welch(y, fs, nperseg=4096, scaling='spectrum', window='hanning')
    
    # Plot the periodograms
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(18, 7))
    fig.suptitle(title, fontsize=18)
    
    ax1.semilogy(fxx, Rxx, label='Espectro de potencia de entrada $x(n)$')
    ax1.set_xlabel('$f$ [Hz]', fontsize=18)
    ax1.set_ylabel('$R_{xx}$', fontsize=18)
    ax1.legend(fontsize=15)
    ax1.grid()
    
    ax2.semilogy(fyy, Ryy, label='Espectro de potencia de salida $y(n)$')
    ax2.set_xlabel('$f$ [Hz]', fontsize=18)
    ax2.set_ylabel('$R_{xx}$', fontsize=18)
    ax2.legend(fontsize=15)
    ax2.grid()
    
    plt.show()

def plot_spectrograms(x, y, fs, title):
    """ Plot the spectrogram of the input signal and the desired response signal
        @param x Input signal
        @param y Desired response signal
        @param fs Sampling frequency
        @param title Title of the plot
    """
    # Size of the windows
    N = 256

    # Plot the spectrogram of the input and the output
    fig, ax = plt.subplots(2, 1, figsize=(20, 15))
    fig.suptitle(title, fontsize=18)
    
    ax[0].set_title('Espectrograma $x(n)$', fontsize=16)
    ax[0].specgram(x, Fs=fs, NFFT=N, window=np.hamming(N))
    ax[0].set_ylabel('$f$ [Hz]', fontsize=16)
    ax[0].set_xlabel('$t$ [s]', fontsize=16)
    ax[0].grid()

    ax[1].set_title('Espectrograma $y(n)$', fontsize=16)
    ax[1].specgram(y, Fs=fs, NFFT=N, window=np.hamming(N))
    ax[1].set_ylabel('$f$ [Hz]', fontsize=16)
    ax[1].set_xlabel('$t$ [s]', fontsize=16)
    ax[1].grid()

    plt.show()

2.1. Estimación espectral

plot_periodograms(x_voice, y_voice, fs, 'Señal de voz')

plot_periodograms(x_music, y_music, fs, 'Señal de música')

plot_periodograms(x_square, y_square, fs, 'Señal rectangular')

plot_periodograms(x_linear_sweep, y_linear_sweep, fs, 'Barrido lineal de frecuencias')

plot_periodograms(x_exponential_sweep, y_exponential_sweep, fs, 'Barrido exponencial de frecuencias')

plot_periodograms(x_noise, y_noise, fs, 'Ruido blanco gaussiano')

2.2. Espectrogramas

plot_spectrograms(x_voice, y_voice, fs, 'Señal de voz')

plot_spectrograms(x_music, y_music, fs, 'Señal de musica')

plot_spectrograms(x_square, y_square, fs, 'Señal rectangular')

plot_spectrograms(x_linear_sweep, y_linear_sweep[:-32000], fs, 'Señal de barrido lineal')

plot_spectrograms(x_exponential_sweep, y_exponential_sweep[:-57000], fs, 'Señal de barrido exponencial')

plot_spectrograms(x_noise, y_noise, fs, 'Señal de ruido blanco')

2.3. Conclusiones

La señal que posea contenido espectral en un mayor rango de frecuencias audibles, será la que más información (en principio) proveerá para estimar la respuesta al impulso del sistema. Esto es porque se sabrá cómo responde el mismo ante estas frecuencias. Por otro lado, es importante además tener en cuenta cómo es la energía en cada frecuencia con la cual se excita al sistema, ya que eso impacta en cómo será la señal deseada sobre la cual se buscará el estimador óptimo. Es decir, a menor energía presente en ciertas bandas, peor será el estimador, aunque exista contenido espectral.

En conclusión, es importante observar que exista contenido en todo el rango audible, pero además que exista una buena distribución de energía. Únicamente viendo los periodogramas, se puede decir que, tanto el barrido lineal de frecuencias como el ruido blanco gaussiano permitirán estimar adecuadamente la respuesta del sistema, ya que poseen contenido espectral para las frecuencias entre 20Hz y 20kHz, con una energía igual para cada modo.

Por otro lado, la señal del barrido lineal de frecuencias permite observar de forma más limpia el comportamiento del sistema a cada frecuencia, mientras que el ruido blanco gaussiano posee un espectro de amplitudes más amplio, y en todo instante el sistema estaría respondiendo a muchas frecuencias. Desde este punto de vista, la mejor estimación se puede hacer con el barrido lineal de frecuencias, ya que se excita al sistema de a un modo a la vez.

3. Elección del modelo

Es necesario determinar el modelo matemático a emplear para la identificación de sistemas. Este modelo será un FIR o modelo de promedios móviles (moving average). En este caso, el modelo asume que el sistema sólo posee ceros, ya que sus polos se encuentran siempre en el origen. Así, el modelo se encuentra completamente parametrizado cuando se definen los coeficientes del mismo. Para esto, es necesario conocer cuántos coeficientes tiene el modelo, es decir, qué orden se utiliza.

Para esto, contando con x(n) e y(n), buscamos la respuesta al impulso $\textbf{h(n)}$ del sistema parlante-habitación-micrófono. El criterio de elección será por inspección del MSE mínimo $J_o$ como función de $M$.

El MSE mínimo se obtiene mediante: \begin{equation} J(\textbf{w}_o) = \sigma_y^2 - \textbf{p}^H\cdot \textbf{w}_o \end{equation} Donde: $\textbf{w}_o = \textbf{R}^{-1}\cdot \textbf{p}$

from scipy import signal

from scipy import linalg

import time

def estimate_correlations(x, y):
    """ Estimates the autocorrelation function of the input signal and the cross correlation
        between the input signal and the desired response of the system.
        @param x Input signal x(n)
        @param y Desired response y(n)
        @return (rxx, rxy)
    """
    # Estimate the autocorrelation of the input signal
    rxx = signal.correlate(x, x, method='fft') / len(x)
    # Estimate the cross correlation between the input and the output signal
    rxy = signal.correlate(y, x, method='fft') / len(x)
    return rxx, rxy

def truncate_correlations(rxx, rxy, lags=None, one_sided=True):
    """ Truncates the autocorrelation and the cross correlation to use only the needed segment.
        @param rxx Estimation of the autocorrelation
        @param rxy Estimation of the cross correlation
        @param lags Amount of lags to each side (positive and negative) to be returned
        @param one_sided Whether it returns only the positive lags or both
        @return (rxx, rxy) Autocorrelation and cross correlation estimations
    """
    # Both estimations have been computed for negative and positive lags. According to what
    # is specified in the function arguments, either positive lags or both will be returned.
    # The following implementation was chosen to be easy to read rather than efficient in
    # amount of code lines.
    if one_sided:
        if lags is not None:
            rxx = rxx[rxx.argmax():rxx.argmax()+lags+1]
            rxy = rxy[len(rxy)//2:len(rxy)//2+lags+1]
        else:
            rxx = rxx[rxx.argmax():]
            rxy = rxy[len(rxy)//2:]
    else:
        if lags is not None:
            rxx = rxx[rxx.argmax()-lags:rxx.argmax()+lags+1]
            rxy = rxy[len(rxy)//2-lags:len(rxy)//2+lags+1]
        
    return rxx, rxy

def optimum_linear_estimator(x, y, order, rxx=None, rxy=None, optimize=True):
    """ Computes the coefficients for the optimum linear estimator using an estimation of
        the autocorrelation and the cross correlation. It also returns an estimation of the
        minimum MSE obtained with the optimum estimator (residual MSE).
        @param x Input signal x(n)
        @param y Desired response y(n)
        @param order Order of the linear estimator
        @param rxx Estimation of the autocorrelation
        @param rxy Estimation of the cross correlation
        @param optimize Whether to use the Levinson-Durbin algorithm or not
        @return (Jo, wo) Residual MSE and optimum coefficients
    """
    # Estimate the autocorrelation and the cross correlation
    if rxx is None or rxy is None:
        rxx, rxy = estimate_correlations(x, y)
    
    if optimize:
        # The solve_toeplitz method (which we'll use since R is toeplitz) takes the first row and column of the matrix R
        # The first row is rxx, the first column is rxx* except for its first element.    
        rxx_col = np.conj(rxx)
        rxx_col[0] = rxx[0]
        wo = linalg.solve_toeplitz((rxx_col, rxx), rxy)
    else:
        # Create the autocorrelation matrix
        R = linalg.toeplitz(rxx)
        wo = np.dot(linalg.inv(R), rxy)

    Jo = y.var() - rxy.dot(wo)
    return Jo, wo

def run_order_analysis(x, y, orders, optimize=False):
    """ Calculates the optimum linear estimator coefficients for each of the given orders
        and estimates the power of the error. Returns the residual MSE and the execution time
        of each estimation.
        @param x Input signal x(n)
        @param y Desired response y(n)
        @param orders List of orders to be computed
        @param optimize Whether to use the Levinson-Durbin algorithm or not
        @return (residuals, execution, coefficients)
    """
    # Compute the estimation of the autocorrelation and the cross correlation
    rxx, rxy = estimate_correlations(x, y)

    # Initialization of numpy arrays and selection of orders
    residuals = np.zeros(len(orders), dtype=np.float32)
    execution = np.zeros(len(orders), dtype=np.float32)
    coefficients = []

    for i, order in enumerate(orders):
        # Improve user experience, reduce anxiety!
        IPython.display.clear_output(wait=True)
        print(f'[{np.round((i+1) / len(orders) * 100, 1)}%] Paso {i+1} de {len(orders)} - Usando M = {order}')
        
        # Starting time of the execution
        start_time = time.time()

        # Compute the estimator
        rxx_truncated, rxy_truncated = truncate_correlations(rxx, rxy, lags=(order-1), one_sided=True)
        Jo, wo = optimum_linear_estimator(x, y, order, rxx=rxx_truncated, rxy=rxy_truncated, optimize=optimize)
        coefficients.append(wo)

        # Save the estimator results
        execution[i] = time.time() - start_time
        residuals[i] = Jo

    print(f'Se estimaron los coeficientes de {len(orders)} órdenes en {np.round(execution.sum(), 2)} segundos')
    return residuals, execution, coefficients

def plot_analysis_results(orders, y, residuals, execution, coefficients, clog=False):
    """ Plots the results of the analysis
        @param orders List with orders evaluated
        @param y Output signal
        @param residuals MSE residuals
        @param execution Execution time of each order
        @param coefficients Coefficients of each order
        @param clog Whether to plot the coefficients in logaritmic scale or not
    """
    # Create the layout for the plots
    fig, ax = plt.subplots(2, 2, figsize=(15, 15))
    ax1 = ax[0][0]
    ax2 = ax[0][1]
    ax3 = ax[1][0]
    ax4 = ax[1][1]
    
    # Plot the minimum MSE obtained for each order
    ax1.plot(orders, [y.var() for order in orders], color='black', label='Peor caso de potencia de salida')
    ax1.plot(orders, residuals, marker='o', color='red', label='MSE residual')
    ax1.set_ylabel('$\hat{J}_o$', fontsize=16)
    ax1.set_xlabel('$M$', fontsize=16)
    ax1.legend(fontsize=15)
    ax1.grid()
    
    # Plot the execution time for each order
    ax2.plot(orders, execution, marker='o', color='blue', label='Tiempo de ejecución')
    ax2.set_ylabel('$t$ [s]', fontsize=16)
    ax2.set_xlabel('$M$', fontsize=16)
    ax2.legend(fontsize=15)
    ax2.grid()

    # Plot the normalized MSE for each order
    ax3.plot(orders, np.round(residuals / y.var() * 100, 2), marker='o', color='green', label='MSE normalizado')
    ax3.set_ylabel('$\\xi$ [%]', fontsize=16)
    ax3.set_xlabel('$M$', fontsize=16)
    ax3.set_yticks(np.linspace(0, 100, 11, dtype=np.int))
    ax3.legend(fontsize=15)
    ax3.grid()
    
    # Plot the coefficients
    if clog:
        ax4.plot(10*np.log10(np.abs(coefficients[residuals.argmin()])), label='Coeficientes')
        ax4.set_ylabel('[dB]', fontsize=16)
    else:
        ax4.stem(coefficients[residuals.argmin()], label='Coeficientes', use_line_collection=True)
    ax4.set_xlabel('$n$', fontsize=16)
    ax4.legend(fontsize=15)
    ax4.grid()

    plt.show()

3.1. Validación del estimador

En esta sección, se realiza una prueba para validar que funciona correctamente el esquema de filtrado lineal óptimo. Para ello, se genera una secuencia de entrada $x(n)$ de ruido blanco gaussiano, y se aplica un filtrado FIR de orden 10. Luego, a la salida se suma otro ruido blanco gaussiano. Finalmente, sobre esta salida $y(n)$ se busca aplicar el esquema de identificación de sistemas.

Los coeficientes del sistema físico a identificar son $h(n) = \frac{1}{10}$ para $n = 0, 1, ..., 9$

# Generate random input signal
x_test = np.random.normal(loc=0, scale=1, size=len(t))

# Obtain the system's response
N = 10
y_test = signal.lfilter(np.ones(N) / N, np.ones(1), x_test)

# Add gaussian white noise with 0,1 variance
y_test += np.random.normal(loc=0, scale=0.01, size=len(y_test))

# Run the analysis
orders = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
test_residuals, test_execution, test_coefficients = run_order_analysis(
    x_test, 
    y_test,
    orders,
    optimize=False
)

[100.0%] Paso 20 de 20 - Usando M = 20
Se estimaron los coeficientes de 20 órdenes en 0.1899999976158142 segundos

# Plot the analysis results
plot_analysis_results(orders, y_test, test_residuals, test_execution, test_coefficients)

Se puede observar que los coeficientes del estimador lineal óptimo se anulan a partir de cierto órden, por ende, el orden óptimo para este modelo sería $M = 10$, lo cual coincide con el orden del sistema físico simulado. Por otro lado, se puede notar que los coeficientes convergen a los del filtro FIR creado. Finalmente, cuando se alcanza el error cuadrátrico medio mínimo, la potencia de error coincide con la potencia del ruido superpuesta a la salida del sistema.

Estas observaciones permiten validar el correcto funcionamiento del estimador desarrollado.

3.2. Análisis de orden para el barrido lineal de frecuencias

Se realiza un análisis del comportamiento del estimador lineal óptimo aplicado a la señal de barrido lineal de frecuencias, utilizando diferentes órdenes para comprobar cuál es el mejor modelo aplicable a este escenario.

# Run the analysis for the linear sweep excitation
orders = [100, 200, 300, 400, 500, 600, 700, 800, 900, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 20000]

residuals, execution, coefficients = run_order_analysis(
    x_linear_sweep[:len(y_linear_sweep[54700:-32000])],
    y_linear_sweep[54700:-32000],
    orders,
    optimize=True
)

[100.0%] Paso 18 de 18 - Usando M = 20000
Se estimaron los coeficientes de 18 órdenes en 2.9100000858306885 segundos

# Plot the analysis results
plot_analysis_results(orders, y_linear_sweep[54700:-32000], residuals, execution, coefficients)

La estimación del error cuadrático medio decae abruptamente a partir de cierto orden, indicando haberse acercado a la condición óptima, y luego para mayores órdenes no hay una mejora apreciable a nivel gráfico. Desde este punto de vista, visualmente se podría considerar necesario modelar al sistema con un orden aproximadamente $M = 10000$.

En la estimación de la respuesta impulsiva $h(n)$, los coeficientes presentan al principio valores nulos o cercanos a nulos. Esto se debe a que las señales de entrada y salida deseada no se encuentran completamente sincronizadas, sino que la salida está retrasada. Este retardo temporal está siendo considerado parte del sistema físico, y por ello nuestra estimación lo modela agregando coeficientes nulos.

Por otro lado, si bien observando la estimación del MSE suponemos que el orden adecuado está alrededor de $M = 10000$. En la respuesta impulsiva no se puede apreciar si efectivamente los coeficientes se anularon para ese orden. Esto podría implicar que en realidad es necesario un orden mayor, aunque la mejoría en el MSE no sea apreciable, el efecto podría ser perceptible a nivel auditivo. Entonces, sabiendo que los ecos en una habitación tienden a provocar un comportamiento decayente exponencialmente.

# Plotting the coefficients in a logarithmic scale
plt.figure(figsize=(10, 5))
plt.plot(10 * np.log10(np.abs(coefficients[residuals.argmin()])), label='Coeficientes')
plt.xlabel('$n$', fontsize=16)
plt.ylabel('dB', fontsize=16)
plt.legend(fontsize=16)
plt.grid()
plt.show()

La caída de la amplitud o envolvente de los coeficientes en escala logarítmica es aproximadamente lineal, como se esperaba, no obstante no se llega a apreciar el punto donde se produce el codo hacia una línea horizontal de amplitudes aproximadamente nulas. Esto sugiere la posibilidad de que existan más coeficientes y que el orden del sistema sea mayor que el considerado. Una solución sería aumentar el orden en el análisis hecho, para ver en qué punto dejan de haber coeficientes no nulos, no obstante con el método utilizado actualmente el costo computacional es muy elevado y requiere muchos recursos de memoria. Entonces, se optó por considerar que el resultado óptimo es con $M = 20000$.

Este resultado deja entrever que el ambiente en donde se realizó las grabaciones es altamente reverberante.

Por otro lado, con respecto al costo computacional del método de cálculo empleado. Se podría reducir el cálculo empleando el algoritmo de Levinson-Durbin, ya que reduce la complejidad computacional necesaria para resolver las ecuaciones normales. De aquí en adelante, se empleará este método.

3.3. Escuchando la respuesta impulsiva del sistema

Se puede escuchar la respuesta impulsiva del sistema parlante-habitación-micrófono.

IPython.display.Audio(coefficients[-1], rate=fs)

Guardamos para más adelante los coeficientes y el MSE

coefs_linear_sweep = coefficients[-1]
mse_linear_sweep = np.round(residuals[-1] / y_linear_sweep[54700:].var() * 100, 2) 

4. Estimación de $h(n)$ para diferentes excitaciones

A continuación se estima $h(n)$ para $M=20000$ con las demás señales presentadas

from numpy.fft import fft

from numpy.fft import fftfreq

def estimate_response(x, y, order, optimize=True):
    """ Calculates the optimum linear estimator coefficients for each given order
        and estimates the error's power. Returns the normalized MSE.

        @param x Input signal x(n)
        @param y Desired response y(n)
        @param order Order to be computed
        @param optimize Whether to use the Levinson-Durbin algorithm or not
        @return (mse, coefficients)
    """
    # Compute the estimation of the autocorrelation and the cross correlation
    rxx, rxy = estimate_correlations(x, y)
    
    # Truncate correlations and estimate wo
    rxx_truncated, rxy_truncated = truncate_correlations(rxx, rxy, lags=(order-1), one_sided=True)
    Jo, wo = optimum_linear_estimator(x, y, order, rxx=rxx_truncated, rxy=rxy_truncated, optimize=optimize)
    
    print(f'Estimated impulse response, MSE is {Jo:.2f}')
    
    # Calculate normalized MSE
    normalized_mse = np.round(Jo / y.var() * 100, 2)
    return normalized_mse, wo

def plot_response(coefficients, fs, estimation, mse, clog=False):
    """ Plots the resulting estimations h(n) and H(exp(jw))
        @param y Output signal
        @param coefs Coefficients
        @param estimation Name for plot title
        @param mse Normalized mse, shown in the title
        @param clog Whether to plot the coefficients in dB scale or not
    """
    # Create the layout for the plots
    fig, ax = plt.subplots(1, 2, figsize=(15, 8))
    ax1 = ax[0]
    ax2 = ax[1]
    fig.suptitle(f'{estimation}. $\\xi$={mse}%', fontsize=18)
    
    # Plot the coefficients
    if clog:
        ax1.plot(10*np.log10(np.abs(coefficients)), label='Coeficientes')
        ax1.set_ylabel('[dB]', fontsize=16)
    else:
        ax1.stem(coefficients, label='Coeficientes', use_line_collection=True)
    ax1.set_xlabel('$n$', fontsize=16)
    ax1.legend(fontsize=15)
    ax1.grid()

    # Plot the FFT
    N = len(coefficients)
    freqs = fftfreq(N, d=1/fs)
    freq_resp = np.abs(fft(coefficients))
    freq_resp_dB = 10*np.log10(freq_resp[:N//2])
    ax2.plot(freqs[:N//2] / 1000, freq_resp_dB, label='$\hat{H}(e^{j\omega})$')
    ax2.set_xlabel('$f\,\,[kHz]$', fontsize=16)
    ax2.set_ylabel('dB', fontsize=16)
    ax2.legend(fontsize=15)
    ax2.grid()
    
    plt.show()

# Chosen order
M = 20000
# Estimations
estimations = ['Voice', 'Music', 'Rectangular', 'Linear Sweep', 'Exponential Sweep', 'White Noise']
estimation_residuals = np.zeros((len(estimations), 1))

4.1 $\hat{h}(n)$ para señal de voz

# Causality and minimum-delay offset
no = 41000

# Calculate the optimal estimator
mse_voice, coefs_voice = estimate_response(
    x_voice[:len(y_voice[no:])],
    y_voice[no:], 
    M, 
    optimize=True
)
estimation_residuals[0] = mse_voice

Estimated impulse response, MSE is 1425930.77

plot_response(coefs_voice, fs, 'Señal de voz', mse_voice)

Como primer conclusión, ha de decirse que no se ha podido estimar la respuesta impulsiva de la habitación, a partir de la excitación de una voz hablando. Los coeficientes mostrados en el gráfico superior dan cuenta de una señal prácticamente compuesta en su totalidad por ruido, y esto se debe a la estocasticidad de la excitación utilizada. El habla humana contiene una gran componente aleatoria, que a la hora de hacer una estimación, dificulta la adaptación del filtro para la identificación de un sistema, aportando una importante componente de ruido a la estimación. Se suma a este problema el hecho de que esta señal no excita a todos los modos en el espectro de la frecuencia, y a los que excita, no lo hace de manera equitativa.

En retrospectiva, habría sido de mayor utilidad emitir un sonido estilo "sshhh" (hacer silencio), el cual se asemejaría más a un ruido blanco (tecnicamente, sería un ruido coloreado), o reproducir un tono con la voz (menos estocástico). De todas formas, ninguna de estas opciones habría solucionado por completo la dificultad del ruido en la estimación, debido a la inherente estocasticidad en la voz humana, ni hubiera dado mejores resultados que otras excitaciones.

4.2 $\hat{h}(n)$ para señal de música

# Causality and minimum-delay offset
no = 43800

# Calculate the optimal estimator
mse_music, coefs_music = estimate_response(
    x_music[:len(y_music[no:])],
    y_music[no:], 
    M, 
    optimize=True
)
# Calculate normalized MSE
estimation_residuals[1] = mse_music

Estimated impulse response, MSE is 105416.90

plot_response(coefs_music, fs, 'Señal de música', mse_music)

A pesar de no ser una excitación idónea, exhibe muy buenos resultados a la hora de estimar la respuesta impulsiva de la habitación. No solo puede apreciarse una forma de respuesta que tiene sentido, sino que logra reducir el MSE normalizado al $ 13.16 \% $. Este comportamiento favorable se debe al buen contenido espectral de esta canción en especial, y sobre todo en frecuencias de mayor interés para el oído humano (ver espectrograma en 2.2), logrando buenos resultados (para el oyente) a la hora de usarse esta respuesta como filtro.

Dentro de sus puntos bajos se destacan los artefactos al inicio y final de la respuesta impulsiva, y el hecho de que la misma no tiende a 0, sino a un nivel constante de ruido. Como con el caso de la voz, hay componente estocástica en esta señal también, y se explica desde allí la presencia de este nivel de ruido, mejorado esta vez por un considerablemente mejor SNR, y el contenido musical de la canción. Además, el periodograma de esta estimación nos muestra un espectro que no coincide del todo con el de la grabación del sweep lineal en la sección 2.1 (considerado la mejor respresentación de la habitación).

4.3 $\hat{h}(n)$ para señal rectangular

# Causality and minimum-delay offset
no = 84000

# Calculate the optimal estimator
mse_square, coefs_square = estimate_response(
    x_square[:len(y_square[no:])],
    y_square[no:], 
    M, 
    optimize=True
)
estimation_residuals[2] = mse_square

Estimated impulse response, MSE is 36479072.23

plot_response(coefs_square, fs, 'Señal rectangular', mse_square)

En este caso no se logra obtener una estimación con sentido de $h(n)$, el MSE normalizado de 96.87% es muy alto y la forma de $h(n)$ no es la esperada.

La razón es que si bien la señal rectangular excita solamente las frecuencias armónicas la fundamental, a la salida hay componentes no despreciables en todas las frecuencias, como se puede observar en el periodograma de la grabación, mostrado en la sección 2.1. Esto es un comportamiento alineal del sistema, causando que el modelo FIR asumido no sea válido.

Hay dos motivos para este comportamiento. En primer lugar, la medición está contaminada con un piso de ruido en todas las frecuencias del espectro, no despreciable al contrastar contra el valor nulo que tienen las frecuencias no armónicas de la fundamental en la excitación, lo cual se aprecia muy claramente en el periodograma de la grabación. Además, gracias a los espectrogramas de las grabaciones de los barridos lineal y exponencial de frecuencia, sabemos que el sistema tiene distorsión armónica, lo cual también afecta en este caso. Sin embargo, el motivo principal es el primero, porque si el segundo fuera muy influyente no se llegaría a una estimación aceptable con ninguna excitación.

4.4 $\hat{h}(n)$ para sweep lineal de frecuencia

Ya fue presentado en la sección anterior, fue el primero en ser estimado. Se muestran con animos de completitud los coeficientes y su Transformada de Fourier.

estimation_residuals[3] = mse_linear_sweep

plot_response(coefs_linear_sweep, fs, 'Sweep lineal de frecuencia', mse_linear_sweep)

Se observa una respuesta en frecuencia muy similar a lo esperado al observar el periodograma de la grabación de esta misma señal. Dado que la excitación tiene la misma energía para todos los modos, la forma del periodograma de la grabación está altamente relacionada con la forma de la respuesta en frecuencia, teniendo en cuenta que la relación entre periodograma y respuesta en frecuencia es el cuadrado del módulo.

Por lo tanto, dado que en la respuesta en frecuencia estimada $\hat{H}(e^{j\omega})$ se observa una forma muy similar a la del periodograma, se concluye que con el barrido lineal de frecuencias se obtuvo una estimación muy buena. Esta valoración es validada por el valor del MSE normalizado, 0.07%, el más bajo de todas las estimaciones.

4.5 $\hat{h}(n)$ para sweep exponencial de frecuencia

# Causality and minimum-delay offset
no = 46200

# Calculate the optimal estimator
mse_exp_sweep, coefs_exp_sweep = estimate_response(
    x_exponential_sweep,
    y_exponential_sweep[no:-57000], 
    M, 
    optimize=True
)
estimation_residuals[4] = mse_exp_sweep

Estimated impulse response, MSE is 41506.50

plot_response(coefs_exp_sweep, fs, 'Sweep exponencial de frecuencia', mse_exp_sweep)

Se observa que la respuesta en frecuencia $\hat{H}(e^{j\omega})$ tiene picos similares a los esperados por la respuesta en frecuencia del sweep lineal, solo que con un desplazamiento en amplitud dependiendo de las bandas de frecuencia, similar a una ecualización. En bajas frecuencias (aproximadamente $ < 4kHz$), la respuesta está amplificada con respecto a las frecuencias medias ($4kHz < f < 15kHz$), mientras que para las altas frecuencias ($>15kHz$), se ve atenuada con respecto a la banda anterior.

Este comportamiento se debe a la distribución de energía de la excitación, observada en su periodograma de la sección 2.1. Como las bajas frecuencias concentran la mayoría de la energía, tienen un mejor SNR y la estimación les da una mayor amplitud que a las frecuencias medias, y mucho mayor que a las altas frecuencias.

Finalmente, se puede observar que la respuesta impulsiva tiene un artefacto en $t=0$ de muy corta duración, el cual es el causante de que para el límite superior de frecuencia $\hat{H}(e^{j\omega})$ sea creciente, al ser un pequeño pico de alta frecuencia.

Esta respuesta también es buena dado que además de asemejarse a la esperada, el MSE normalizado es bastante bajo, de 0.29%.

4.6 $\hat{h}(n)$ para ruido blanco gaussiano

# Causality and minimum-delay offset
no = 38000

# Calculate the optimal estimator
mse_noise, coefs_noise = estimate_response(
    x_noise[:len(y_noise[no:])],
    y_noise[no:], 
    M, 
    optimize=True
)
estimation_residuals[5] = mse_noise

plot_response(coefs_noise, fs, 'Ruido blanco gaussiano', mse_noise)

Estimated impulse response, MSE is 80741.02

Finalmente, se evalúan los resultados obtenidos al excitar con ruido blanco. A priori, siendo una señal que excitaba de forma pareja a todos los modos, se esperaba poder lograr una buena estimación de la respuesta impulsiva. Sin embargo, se observa en el gráfico un importante nivel de ruido presente en los coeficientes. Al igual que con el caso de la voz humana, la estocasticidad de la señal va en detrimento del funcionamiento del filtro adaptativo, y justifica este alto nivel de ruido. De todas formas, se logran mejores resultados que en el caso de la voz, dado que no solo el MSE pudo ser reducido en un pequeño porcentaje, sino que se llega a apreciar el inicio de la respuesta impulsiva por encima del ruido.

4.7 Conclusión

Los mejores resultados fueron obtenidos con los sweep de frecuencias, siendo ligeramente superior el lineal por su distribución equitativa de la energía. Estas señales con un comportamiento más predecible y menos aleatorio, lograron respuestas impulsivas con muy bajos niveles de ruido, lo cual no fue el caso de otras señales con mayor componente estocástica, como el ruido blanco, la voz y la señal de música. La cuadrada presentó la performance más baja, aunque esto era esperado por el hecho de no excitar a todos los modos.

5. Comparación con la predicción

En la Sección 2 se predijo que la mejor estimación iba a venir dada por el sweep lineal de frecuencias, lo cual terminó sucediendo, en términos matemáticos. Además se vio en los resultados que el sweep exponencial de frecuencia, si bien no cumple la condición pedida de que la energía sea igual para todas las frecuencias, logró una estimación bastante buena.

Contrario a lo esperado por el hecho de tener contenido espectral en todas las frecuencias y de forma pareja, se vió que el ruido blanco produce una mala estimación, lo cual se debe a su naturaleza estocástica, como se explica en la sección anterior.

6. Elección de mejor estimación y filtrado

6.1 Elección de mejor estimación

Junto con lo comentado en las secciones anteriores, se muestra a continuación una comparación de los MSE de cada estimación para elegir la mejor $\hat{h}(n)$

# Plot MSE for each estimation
plt.scatter(estimations, estimation_residuals)
plt.xticks(rotation=90, fontsize=12)
plt.xlabel('', fontsize=15)
plt.ylabel('$\\xi$ [%]', fontsize=14)
plt.title('MSE para cada estimación', fontsize=16)
plt.grid()
plt.show()

Los mejores MSE son obtenidos con los sweeps de frecuencia, siendo ligeramente mejor entre ambos el sweep lineal. Esto se corresponde con lo observado y explicado en el dominio de la frecuencia, ya que el sweep lineal estima correctamente hasta aproximadamente $17.5kHz$, mientras que el sweep exponencial se limita a $10kHz$. Por lo tanto, un mayor rango de frecuencias estimadas correctamente se corresponde con un mejor MSE.

En consecuencia, la estimación con el sweep lineal de frecuencias es considerada la mejor numéricamente.

6.2 Evaluación perceptual de la señal filtrada

A continuación se filtra la señal de música tanto con la mejor estimación como con la estimación del sweep exponencial de frecuencias.

# Filtramos con el sweep lineal
y_music_filt_l = signal.lfilter(coefs_linear_sweep, np.ones(1), x_music)

# Escuchamos con el sweep lineal
IPython.display.Audio(y_music_filt_l, rate=fs_y_linear_sweep)

Al filtrar la señal de música con la estimación obtenida por el sweep lineal se percibe como si hubiera sido grabado en una habitación similar a la utilizada. Sin embargo, la distribución energética en frecuencias, genera una percepción más aguda de la señal, es decir con frecuencias mas altas. Por lo tanto, se prueba filtrar ahora el audio con la exitación sweep exponencial que también tenía un MSE muy bajo para ver los resultados.

# Filtramos con el sweep exponencial
y_music_filt_e = signal.lfilter(coefs_exp_sweep, np.ones(1), x_music)

# Escuchamos con el sweep exponencial
IPython.display.Audio(y_music_filt_e, rate=fs_y_exponential_sweep)

6.3 Conclusiones finales

Si bien se mencionó en la sección 6.1 que la señal que mejor estima la respuesta impulsiva es la del sweep lineal, debido a los resultados numéricos obtenidos, si se procede a escuchar y comparar los resultados sobre el audio de la música, se logra apreciar una notable diferencia entre lo estimado mediante el sweep lineal y el exponencial. El audio del sweep exponencial es, a diferencia de lo esperado previamente, más similar al audio original. Esto se explica desde el conocimiento de que el oído humano tiene un comportamiento logarítmico, o en otras palabras, es más sensible a las frecuencias donde la relación señal a ruido del sweep exponencial es mas alta. Puede observarse en los espectrogramas de la sección 2.2, que la banda de frecuencia hasta los $ 10 kHz $ posee un claro mayor SNR que la banda de $ 10 kHz $ a $ 20 kHz $, logrando así mejores resultados auditivos que con el sweep lineal.